%get the files you need -- 40 to 52 weeks (10 - 12 months)
mainpath = '/Volumes/Looxcie/1C_LooxcieAdult/03_JPG_AllMasters';
allfolders = dir('/Volumes/Looxcie/1C_LooxcieAdult/03_JPG_AllMasters/**/*.jpg');
%% sample every 5 s (every 25 images)
imnum = zeros(size(allfolders,1),1);
keep = zeros(size(allfolders,1),1);
for i = 1:size(allfolders,1)
    bits = split(allfolders(i).name, "_");
    num = bits{2};
    num2 = num(1:6);
    num3 = str2num(num2);
    imnum(i) = num3;
    if mod(num3,25) == 0
        keep(i) = 1;
    end
end
framelist = allfolders(keep == 1,:);
framelist = struct2table(framelist);
framelist = unique(framelist,'rows');
writetable(framelist, '/Volumes/Looxcie/1C_LooxcieAdult/adultframelist2.csv');
%% set processing paths, etc.
framelist = readtable('/Volumes/Looxcie/1C_LooxcieAdult/adultframelist2.csv');
sfPath = '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/contrast and sf/infant clutter code';
pyrPath = '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/orientation';

% low frequency cut-off, removes frequencies below this point in amplitude
% plot and fit
cutoffLF = 0.1; 
FreqPlot = 30;   % Number of freqs to plot in quick first view of amp spectrum (fig 3)
% cpi to cpd division constant
cpdCon = 41;  % vertical image size in degrees of Looxcie (480 pixels)

% get all files in selected folder, recursively
numImages = size(framelist,1);
 
% Pick one of these 3 window forms to apply to images
windowType = 'cos';
%% set up params
% read in first image to get dimensions of all of them
fn = fullfile(framelist.folder{1}, framelist.name{1});
currImage = imread(fn);

% dimension first to make sure there is enough memory
imageDims = size(currImage);
imageDims = imageDims(1:2);

% struct to store parameters and data
params.impath = framelist.folder{:};
params.fileNames = repmat("",numImages,1);
params.numImages = numImages;
params.imageDims = imageDims;

% store fourier images and vectors
params.radialFreqMag = zeros(numImages,min(imageDims)/2);
params.rmsImageCon = zeros(numImages,1);
params.edges_full = zeros(numImages,1);
params.edges_center = zeros(numImages,1);
params.totalPower = zeros(numImages,1);
params.domOr = zeros(numImages,1);
%% run images through for loop
for fileCounter = 1:numImages
    fn = fullfile(framelist.folder{fileCounter}, framelist.name{fileCounter});
    currImage = imread(fn);
        
    imageDims = size(currImage);
    imageDims = imageDims(1:2);
    
    h = imageDims(1)/2;
    w = imageDims(2)/2;
    
    I2 = rgb2gray(double(currImage)./256);
    
    edgez1 = edge(I2,'canny',[0.10 0.25]);
    params.edges_full(fileCounter) = sum(edgez1, 'all');
    
    edgez1(1:5,:) = 0;
    edgez1((h-4):h,:) = 0;
    edgez1(:,1:5) = 0;
    edgez1(:,(h-4):h) = 0;
    
    params.edges_in(fileCounter) = sum(edgez1, 'all');
    
    cd(sfPath);
    addpath(genpath('/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/contrast and sf/infant clutter code/MATLAB'))

    
    if strcmp(windowType,'cos')
        windowWidth = round(min(imageDims)*0.8);
        imageWindow = cosine2D(min(imageDims),windowWidth);
        %need tukeywin - from signal processing toolbox
    elseif strcmp(windowType,'gauss')
        numGaussSd = 6;                         % number of sd of gaussian window
        gaussSd = min(imageDims)/numGaussSd;    % sd of window, in pixels
        imageWindow = dogauss2d(min(imageDims),gaussSd,0,gaussSd);
    else
        imageWindow = ones(min(imageDims));
    end

    % get the central square region of the gray image for fourier analysis
    Looxcie = mean(currImage,3);
    Looxcie = getcenter(Looxcie,min(imageDims));
    % window image (USE THIS ONE: Window in contrast rather than luminance)
    LooxcieMeanLum = mean(mean(Looxcie));
    Looxcie = ((Looxcie-LooxcieMeanLum).*imageWindow)+LooxcieMeanLum;
    % window image in luminance

    currFftImage = fft2(Looxcie);                   % take the fft
    currMagImage = fftshift(abs(currFftImage));
    currMagImage = currMagImage./min(imageDims)^2;      % divide by the number of pixels to scale it correctly

    % generate the summary amplitude spectrum
    ftmatrix =currMagImage;
    % make radii coordinates
    c = size(ftmatrix,1)/2+1;
    rmat = (1:size(ftmatrix,1))'*ones(1,size(ftmatrix,1))-c;
    cmat = rmat';
    radiusmat = round(sqrt(rmat.^2 + cmat.^2));
    freq = 0:max(size(radiusmat))/2-1;
    freqamp = zeros(size(freq));
    for i = 1:length(freq)
        currfreq = freq(i);
        currlocs = intersect(find(radiusmat>=currfreq),find(radiusmat<currfreq+1));
        freqamp(i) = sum(ftmatrix(currlocs(:)));  %% sum all of the energy at locations with that frequency
    end

    % output:
    params.freq = freq;
    % convert the frequencies from units of cycles per image (cpi) to cycles per degree (cpd)
    params.freq = params.freq./cpdCon;

    params.radialFreqMag(fileCounter,:) = freqamp;
    params.rmsImageCon(fileCounter) = sqrt(mean((Looxcie(:)-mean(Looxcie(:))).^2));
    % compute total power in fourier spectrum 
    params.totalPower(fileCounter) = sum(params.radialFreqMag(fileCounter,:));

    % calculate the slope of the function relating frequency to amplitude
    % (in log-log space)
    keeplocs = find(params.freq>=cutoffLF); % locations of frequencies to keep
    inputVals = [params.freq(keeplocs)',params.radialFreqMag(fileCounter,keeplocs)']; % x and y values for fit
    showFit = 0; % toggles plotting the data (1 = show, 0 = suppress)
    initialParams = [max(params.radialFreqMag(keeplocs)),1.5]; % initial parameter estimates
    [params.AmpSpecSlope(fileCounter,1:2), ~ ] =  invffitfn(inputVals,showFit,[],initialParams); 
    
    im = rgb2gray(currImage);
    nOri = 16; % number of orientations to perform the calculations on.

    % Some characteristics of the image.
    nPixels = numel(im);
    dims = size(im); % size of the image
    ctr = ceil((dims+0.5)/2); % center of the image

    % Calculate Fourier Transform and corresponding frequency axis.
    Fs = size(im,1);
    ff = -Fs/2:(Fs/2-1);
    imdft = abs(fftshift(fft2(im)))/nPixels; % amplitude (normalized by nSamples)

    % pre-calculate meshes for computing the filters
    [XX,YY] = meshgrid( ([1:dims(2)]-ctr(2))./(dims(2)/2), ...
        ([1:dims(1)]-ctr(1))./(dims(1)/2) );
    TH = atan2(YY,XX);
    rad = sqrt(XX.^2 + YY.^2);
    rad(ctr(1),ctr(2)) =  rad(ctr(1),ctr(2)-1);  % this line is getting rid of zero
    log_rad  = log2(rad);

    % pre-calculate position for the spatial frequency filter generation
    sf_positions = -[1/16,.5:.5:(round(log2(min(dims)))/2 +.5)]; % top/bottom edges of the filters
    sf_edges = size(im,1)*2.^(sf_positions-1);
    nSF = length(sf_positions)-1;

    % which orientations are we calculating this at.
    orientation = 90:-(180/nOri):-78.75; % Erin was right about the filter order and their orientations.

    % pre-allocate some storage variables for the loop.
    ssq = zeros(size(im));
    cnt = 0;
    amp = zeros(nSF,nOri);
    mid_freq = zeros(nSF,1);
    
    restoredefaultpath
    cd(pyrPath);
    addpath('matlabPyrTools-master');

    % the loop
for ee = 1:nSF
    
        % RADIAL FILTER
        
        % upper bound
        if sf_positions(ee) == 1/16, width = 1/16;
        else, width = 1/4; 
        end
        position = sf_positions(ee);
        upper_radial_mask = generateRadialFilter(width,position,log_rad,true);

        % lower bound
        if ee < length(sf_positions)
            width=1/4;
            position = sf_positions(ee+1);
            lower_radial_mask = generateRadialFilter(width,position,log_rad,false);
        else
            lower_radial_mask = ones(size(upper_radial_mask)); 
        end
        
    radial_mask = upper_radial_mask.*lower_radial_mask;
    for aa=1:nOri

        % ANGLE FILTER
        anglemask = generateAngleFilter(aa,nOri,TH);
        mask = abs(anglemask).*radial_mask;

        cnt=cnt+1;
        
        % this can be used to demonstrate that the sum of squares of all
        % filters adds to 1 (except for at DC and high-frequency leftovers)
        ssq = ssq+mask.^2; 
        
        % selecting the amplitude associated with this filter.
        amp(nSF-ee+1,aa) = sum((mask(:).*imdft(:)));

    end
    
    % calculating the mid-point frequency associated with this SF band
    mid_freq(ee) = mean(sf_edges(ee:ee+1));
end

    amps_by_or = sum(amp);
    tot = sum(amps_by_or);
    amp_prop = amps_by_or/tot;
    A = sort(amp_prop,'descend');

    params.domOr(fileCounter) = A(1);
    
    if mod(fileCounter,100) == 0
               disp(fileCounter)
               save('/Volumes/Looxcie/1C_LooxcieAdult/adult_params.mat','params')
    end
    if fileCounter == numImages
               disp(fileCounter)
               save('/Volumes/Looxcie/1C_LooxcieAdult/adult_params.mat','params')
    end
end
%% output csv file
newDF = array2table(params.radialFreqMag);
cnames = strings(length(params.freq),1);
for i = 1:length(params.freq)
    cnames(i)=sprintf('%d',params.freq(i));
end

newDF.Properties.VariableNames = cnames;
newDF.fileNames = params.fileNames;

%% extract sf - output csv
params.edges_in2 = params.edges_in.';
low = NaN(length(params.domOr),1);
mid = NaN(length(params.domOr),1);
high = NaN(length(params.domOr),1);
for i = 1:length(params.domOr)
   L = sum(params.radialFreqMag(i,1:80));
   M = sum(params.radialFreqMag(i,81:160));
   H = sum(params.radialFreqMag(i,161:240));
   tot = sum(params.radialFreqMag(i,1:240));
   low(i) = L/tot;
   mid(i) = M/tot;
   high(i) = H/tot; 
end
% add ampslope to this
SF_sum = array2table([params.rmsImageCon(1:length(params.domOr)),params.edges_in2,params.edges_full(1:length(params.domOr)),params.AmpSpecSlope(:,2)]);
SF_sum.Properties.VariableNames = {'rmsImageCon','edges_in','edges_full','AmpSlope'};% add ampslope to this
SF_sum.domOr = params.domOr;
SF_sum.filename = framelist.name(1:length(params.domOr));
SF_sum.low = low;
SF_sum.mid = mid;
SF_sum.high = high;

writetable(SF_sum, '/Volumes/Looxcie/1A_LooxcieBloomington/05_coding/13_VisualStatistics_2022/contrast and sf/properties_adults.csv')
